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Abstract. I describe ongoing work developing Bayesian methods for flexible modeling of arrival 
time series data without binning, aiming to improve detection and measurement of X-ray and 
gamma-ray pulsars, and of pulses in gamma-ray bursts. The methods use parametric and semi- 
parametric Poisson point process models for the event rate, and by design have close connections 
to conventional frequentist methods currently used in time-domain astronomy. 
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Measuring the arrival times, directions, and energies of individual quanta — photons or 
l-H ■ particles — potentially provides the finest possible resolution of dynamical astronomical 
phenomena, particularly for high-energy sources producing low detectable fluxes. The 
O |. simplest methods for signal detection and measurement bin the data for statistical or 
computational convenience (e.g., to allow use of asymptotic Gaussian approximations or 
to enable fast Fourier decomposition with an FFT). But methods that instead directly 
analyze the event data without binning can detect weaker signals and probe shorter time 
C3 , scales than methods that require binning. 

The ongoing work I briefly describe here is motivated by studies of X-ray and gamma- 
ray pulsars, producing periodic signals, and gamma-ray bursts, producing chaotic signals 
typically comprised of multiple overlapping pulses. In the former case there may be less 
than one event per period (particularly in energy- resolved studies); in the latter, time 
scales as short as milliseconds are relevant, and detected photons are sparse at high 
energies. Both phenomena motivate development of data analysis techniques that can 
milk every hard-won event for what it is worth. 

For simplicity we here focus just on the arrival time data (also known as time-tagged 
' event (TTE) data), presuming the events being analyzed have been selected to have 
t-H . directions consistent with an origin from a single source, and that energy dependence of 
any putative signal is not significant (so the signal's temporal signature is not corrupted 
• rH , by ignoring event energies). We can represent the data as points on a timeline, as in Fig.[TJ 
The dots denote events at times detected within small time intervals, 5t, representing 
the instrumental time resolution. The empty intervals, denoted Atj, are informative; 
seeing no events in an observed interval provides a constraint on the signal, in contrast 
to simply not observing during the interval. 

Conventional approaches to detecting signals in such data (binned or unbinned) adopt 
a frequentist hypothesis testing approach: one devises a test statistic that measures de- 
parture of the data from the predictions of an uninteresting "null" model, and uses it to 
see if the null model may be safely rejected (implying an interesting signal is present). 
No explicit signal model is needed to define such a test. 
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Figure 1. Arrival time series depicted as points on a timeline, with detection and 
nondetection intervals noted. 

Here I describe methods developed using the Bayesian approach, where one compares 
the null model to explicit alternative models describing interesting signals. One of my 
motivations is to show how conventional "alternative free" test statistics arise in this 
framework (exactly or approximately) . This illuminates implicit assumptions underlying 
conventional methods; more constructively, it provides a framework where generalizations 
of the implicit models may lead to new methods. There are several other virtues of 
adopting a Bayesian approach to these problems, both pragmatic and conceptual; I briefly 
note a few below. A more extensive but still introductory discussion is in Loredo (2011). 

We will model the data with a nonhomogeneous Poisson point process in time. A model, 
M, specifies an intensity function (event rate) r(t; V) that depends on the model's param- 
eters, V . Parametric models have a parameter space of fixed dimension; e.g., a periodic 
signal model will typically have frequency, amplitude, and phase parameters, and possi- 
bly additional parameters describing the light curve shape. Nonparametric models have a 
parameter space whose (effective) dimension may grow with sample size, adapting to the 
data; it may formally be infinite-dimensional. Semiparametric models have a parameter 
space with a fixed-dimension part (e.g, the frequency and phase of a periodic model) and 
a nonparametric part (e.g., an adaptive light curve shape). 

The data drive Bayesian inferences via the likelihood function, the probability for the 
data, D = ({U}, {Atj}), given values for the parameters. Referring to Fig. [I] we build the 
likelihood function by calculating the product of Poisson counting probabilities for zero 
counts in the empty intervals, and one count in each detection interval. The zero-count 
probabilities are of the form exp[— J A dtr(t)\, and the one-count probabilities are of the 
form [r(ti)5t] exp[— r(t)5t] (presuming 5t is small so that r(t)8t <C 1). Thus the likelihood 
function is 

N 

Y[r(U)St, (0.1) 
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where T denotes the full observing interval and N is the number of detected events. To 
go further, we must specify specific rate models and priors for the model parameters. To 
fit a particular model, we use Bayes's theorem to calculate a posterior probability for the 
parameters, p(V\D,M) = p(V\M)Cm('P)/Zm> where Zm is a normalization constant, 
given by the integral of the product of the prior probability density, p(V\M), and the 
likelihood function. To detect a signal, we instead use Bayes's theorem on a hypothesis 
space including one or more models for interesting signals, and the null model (here, 
a constant rate model with a single parameter, the amplitude, A, with r(t) = A). In 
this space, the normalization constant for a particular model, Zm, plays the role of the 
likelihood for the model (as a whole); Zm is often called the marginal likelihood for the 
model, where "marginal" refers to the integration over V used to calculate its value. 

As a simple starting point, consider a model where the logarithm of the rate is pro- 
portional to a sinusoid plus a constant; the logarithm guarantees that the rate itself is 
non-negative. Then we may write the rate as r(t) = Aexp[KCOs(w£ — 4>)]/Io(n), where A 
is the time- averaged rate, and Iq(-) denotes the modified Bessel function of order (this 
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normalizes the exponential factor so that A is the time-average). Light curves with this 
shape have a single peak per period; its width (equivalently, the duty cycle) is determined 
by the concentration parameter, k, with large values corresponding to sharp peaks, and 
k = corresponding to a constant rate. To estimate the frequency and concentration, 
we calculate the posterior for all four parameters, p(A,u},K,<j>\D,M), and marginalize 
(integrate) over A and <j>. Adopting a flat prior for the phase, and nearly any prior for 
A that is independent of the other parameters, we find a marginal posterior probability 
density for frequency and concentration proportional to Iq[kR(lo)}/[Iq(k)\ n , where R(uj) 
is the Rayleigh statistic, given by 



Estimation of to alone, accounting for uncertainty in all other parameters, is found by 
further integrating over k; this is easy to do numerically. Detection requires calculating 
the marginal likelihood, corresponding to a further integration over w. This must be done 
numerically. It can be time-consuming for blind searches, but not significantly more so 
than the kind of frequency grid searching done with conventional tests. 

The Rayleigh statistic was invented for the well-known (frequentist) Rayleigh test 
for periodic signals in arrival time series (see Lewis 1994 for a review of the Rayleigh 
test and other frequentist period detection methods mentioned below). The quantity 
2R 2 {uj) is called the Rayleigh power; it is the point process analog of the periodogram 
or Fourier power spectral density. From a Bayesian point of view, the Rayleigh test 
implicitly assumes periodic signals may be well-modeled by log-sinusoid rate functions. 
Notably, there is no parameter corresponding to k in the Rayleigh test; also, in practice 
it is known to work well only for smooth light curves with a single broad peak per period. 
Such light curves correspond to values of k near unity, another implicit assumption of the 
Rayleigh test. These results indicate that one can implement Bayesian period searches 
using conventional computational tools already at hand for the Rayleigh test. They also 
indicate that explicit consideration of the k parameter may lead to procedures more 
sensitive to sharply-peaked light curves than the conventional Rayleigh. 

Many pulsars have light curves with two or more peaks per period. This suggests 
generalizing the log-sinusoid model to a log-Fourier model, with the logarithm of the 
rate proportional to a sum of harmonic sinusoids. Adopting a finite sum of m harmonics 
with concentration parameters n k [k — 1 to to, with the fundamental corresponding 
to k = 1), we may proceed analogously to the above analysis. The larger number of 
phases and concentration parameters thwarts analytical integration; in an approximate 
treatment the posterior distribution for frequency and concentration is proportional to 
exp[5(w)] with S(uj) = X^feLi KkR(ku). The frequentist test generalizing the Rayleigh 
test to multiple harmonics is the test, with Z^co) — 2 Y^k=i R 2 {kuj), a sum of 
Rayleigh powers at harmonics. Notably, the Bayesian analysis uses the sum of R(koo) 
values ("harmonic Rayleigh amplitudes") rather than powers. Note that, for = 1, 
S 2 (uj) = Z^ n {uj)/2 + J2k ~52j^k R(jw)R(kuj); that is, S 2 contains information not in Z 2 ^. 
Roughly speaking, Z 2 ^ corresponds to incoherently summing power in harmonics, but the 
quantity arising in the Bayesian treatment of a harmonic model instead sums amplitudes, 
accounting for phase information ignored by Z 2 ^. 

A popular frequentist period detection method that aims to be sensitive to periodic 
signals of complex shape is x 2 epoch folding (x 2 -EF) . One folds the arrival times modulo 
a trial period to produce a phase, 6i, for each event, with 6i in the interval [0, 2n]. For 
a constant signal, the phases should be uniformly distributed. The x 2 -EF method bins 
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the phases into B bins and uses Pearson's x 2 to test consistency of the binned phases 
with a uniform distribution. This motivates a Bayesian model with a piecewise-constant 
rate function with B steps per period. Gregory & Loredo (1992; GL92) analyzed this 
model, with the rate in a bin given by Afk, where A is the average rate, and the shape 
parameters, fk (k = 1 to B), specify the fraction of the rate attributed to each bin, 
with ~^2 k /fc = 1. GL92 assigned a constant prior to the shape parameters, based on 
the intuition that this spreads probability across all possible shapes. After marginalizing 
over the shape parameters, the posterior for frequency and phase is inversely proportional 
to the multiplicity of the set of counts of events in phase bins. In a large-count limit, 
this is approximately exp(— x 2 /2), providing a tie to % 2 -EF. The method performed 
impressively, detecting an X-ray pulsar where the Rayleigh test failed, and performing 
well in a simulation study by Rots (1993) comparing it to other methods. 

Despite these successes, there is room for improvement in the GL92 analysis, for a sur- 
prising reason. The constant prior adopted in GL92 in fact does not spread probability 
over all possible shapes. As the number of bins increases, the constant prior puts ever 
larger probability in the neighborhood of flat models, making it harder than necessary 
to detect narrow peaks. The reason is a "curse of dimensionality" known as concen- 
tration of measure: a multidimensional distribution built out of the product of broad 
one-dimensional distributions with finite moments concentrates its probability in a de- 
creasing volume of parameter space as dimension increases. Concentration can be avoided 
by letting the parameters of the one-dimensional component distribution vary with the 
target dimension. A theoretically appealing way to do this is to require divisibility of the 
prior, e.g., the four-bin prior should reduce to the two-bin prior if we create a two-bin 
model out of the combination of bins 1 and 2, and 3 and 4. A divisible Dirichlet distri- 
bution prior accomplishes this, and improves sensitivity to sharply peaked light curves 
so long as any constant background component is small (Loredo 2011). 

Extending the construction to functions described with an infinite number of bins 
or points leads one to consider infinitely divisible priors for nonparametric functions: 
Gaussian process priors for curve fitting, Dirichlet process priors for modeling probabil- 
ity densities, and Levy process priors for modeling Poisson intensities. With a team of 
statistician and astronomer colleagues, I am developing methods using priors built on 
Levy processes for modeling pulses in gamma-ray bursts. This approach can quantify 
uncertainty even in the regime where pulses are highly overlapping. Its implementation 
involves compound Poisson processes, as arise in simple models of accumulation of rain, 
where drops with a distribution of sizes fall radomly over a region of space. And so 
Bayesian modeling of arrival time series has led us from sines to steps to droplets. 
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